################################################################################
#target_group_bkgd
source('/home/jc148322/scripts/libraries/cool_functions.r')
library(SDMTools)
wd="/home/jc148322/Tammy/Target_bkgd/"; dir.create(wd); setwd(wd) #define and set the working directory
climate.dir='/home/jc165798/Climate/CIAS/AWT/250m/baseline.76to05/bioclim/'
occur.file = "/home/jc149336/Maxent/Layers/Occur/WTBIOCLIM_Occurrence_noveg.csv" #define the file with the species observations
maxent.jar = "/home/jc149336/Maxent/maxent.jar" #define the location of the maxent.jar file
species = c('Chaejobe','Chalgoul','Chalnigr','Hippater', 'Hippdiad', 'Hippsemo', 'Macrgiga',
'Miniaust', 'Minimoo', 'Mormbecc', 'Murflor', 'Myotmacr', 'Nyctbifa', 'Nyctgoul',
'Phonpapu', 'Pipiadam', 'Rhinmega', 'Rhinrobe', 'Saccflav', 'Saccsacc', 'Scotgrey',
'Scotorio', 'Scotruep','Scotsanb', 'Taphaust', 'Taphtrou', 'Vesppumi', 'Vesptrou') #define the species of interest

###prepare the data
occur = read.csv(occur.file,as.is=TRUE) #read in the occur file
clim.vars = grep('bioclim',colnames(occur),value=T) #define the climate variables of interest
occur = occur[,1:3]
occur = unique(na.omit(occur)) #remove any missing data and keep on unique data

for (clim.var in clim.vars) { cat(clim.var,'\n') #cycle through each of hte climate variables
	occur[,clim.var] = extract.data(cbind(occur$long,occur$lat), read.asc.gz(paste(climate.dir,clim.var,'.asc.gz',sep=''))) #append the climate data
		
}
bkgd=occur; bkgd$spp='bkgd'
write.csv(bkgd,'bkgd.csv',row.names=FALSE) #write out your target group background


###cycle through each species and submit jobs to be modelled
for (spp in species) { cat(spp,'\n')
	toccur = occur[occur$spp==spp,] #get the observations for the species
	spp.dir = paste(wd,'models/',spp,'/',sep='') #define the species directory
	dir.create(paste(spp.dir,'output',sep=''),recursive=TRUE) #create the species directory
	write.csv(toccur,paste(spp.dir,'occur.csv',sep=''),row.names=FALSE) #write out the file
	 
	#create the shell script which will run the modelling jobs
	zz = file(paste(spp.dir,'01.',spp,'.model.sh',sep=''),'w') #create the shell script to run the maxent model
		cat('#!/bin/bash\n',file=zz)
		cat('cd ',spp.dir,'\n',sep='',file=zz)
		cat('source /etc/profile.d/modules.sh\n',file=zz) 
		cat('module load java\n',file=zz)
		cat('java -mx2048m -jar ',maxent.jar,' -e ',wd,'bkgd.csv -s occur.csv -o output nothreshold nowarnings novisible replicates=10 nooutputgrids -r -a \n',sep="",file=zz) #run maxent bootstrapped to get robust model statistics
		cat('cp -af output/maxentResults.csv output/maxentResults.crossvalide.csv\n',file=zz) #copy the maxent results file so that it is not overwritten
		cat('java -mx2048m -jar ',maxent.jar,' -e ',wd,'bkgd.csv -s occur.csv -o output nothreshold outputgrids plots nowarnings  responsecurves jackknife novisible nowriteclampgrid nowritemess writeplotdata -P -J -r -a \n',sep="",file=zz) #run a full model to get the best parameterized model for projecting
	close(zz)
	setwd(spp.dir); system(paste('qsub -m n 01.',spp,'.model.sh',sep='')); setwd(wd) #submit the script
}
################################################################################
###To project the models###
##Wherever you see output/ascii/  I would change that to output/ascii_aus/  or ascii_png##
##NEW
work.dir='/home/jc148322/Tammy/Target_bkgd/models/'
mxe.dir = '/home/jc165798/Climate/CIAS/AWT/250m/bioclim_mxe/'
maxent.jar = "/home/jc149336/Maxent/maxent.jar"
proj.list = list.files(mxe.dir) #list the projections
proj.list = proj.list[c(grep('current',proj.list),intersect(grep('RCP',proj.list),c(grep('2035',proj.list),grep('2055',proj.list),grep('2085',proj.list))))] #subset it to simply current and 2080 data

species = list.files(work.dir) #get a list of species

#cycle through each of the species
for (spp in species) { #spp=species[1]
	spp.dir = paste(work.dir,spp,'/',sep=''); setwd(spp.dir) #set the working directory to the species directory

	zz = file(paste('02.',spp,'project.models.sh',sep=''),'w') ##create the sh file
		cat('#!/bin/bash\n',file=zz)
		cat('cd ',spp.dir,'\n',sep='',file=zz)
		cat('source /etc/profile.d/modules.sh\n',file=zz)
		cat('module load java\n',file=zz)
		dir.create('output/ascii_WTcurrent/',recursive=TRUE) #create the output directory for all maps
		#cycle through the projections
		for (tproj in proj.list) cat('java -mx1024m -cp ',maxent.jar,' density.Project ',spp.dir,'output/',spp,'.lambdas ',mxe.dir,tproj,' ',spp.dir,'output/ascii_WTcurrent/',tproj,'.asc fadebyclamping nowriteclampgrid\n',sep="",file=zz)
		cat('gzip ',spp.dir,'output/ascii_WTcurrent/*asc\n',sep='',file=zz)#Try not to zip
	close(zz)

	#submit the script
	system(paste('qsub -m n 02.',spp,'project.models.sh',sep=''))
	}
